##################################################################
###                    Make map of incidents                   ###
##################################################################

## Use file "crosssectional.dta"

library(foreign)
library(plyr)
library(rgeos)
library(maptools)
library(gpclib)  
library(ggplot2)
library(Cairo)
library(ggmap)
library(scales)
library(RColorBrewer)
library(rgdal)

## Don't use original file
myvars <- c("FIPS", "Total")
incident_map <- crosssectional[myvars]
 
## Load map shapefile from local working directory
## First, unzip cb_2015_us_county_5m.zip in local working directory
library(raster)
countymap <- shapefile("cb_2015_us_county_5m.shp")

## Adjust FIPS as unit, then recreate after fortifying
require(sp)
countymap$FIPS <- as.numeric(as.character(countymap$GEOID))
countymap <- fortify(countymap, region = "FIPS")
countymap$FIPS <- as.numeric(countymap$id)

## Merge with incident count
incident_map$FIPS <- as.numeric(as.character(incident_map$FIPS))
graph_incidentmap <- merge(countymap,incident_map,by="FIPS", all=F)

## Get map data for US  states
state_map <- map_data("state")
west_coast <- subset(state_map, region %in% c("california", "oregon", "washington","nevada",
                                              "new mexico","colorado","utah","arizona","montana",
                                              "idaho","wyoming"))
county_map <- map_data("county")
cborders <- subset(county_map, region %in% c("california", "oregon", "washington","nevada",
                                             "new mexico","colorado","utah","arizona","montana",
                                             "idaho","wyoming"))

## Need to have zeros become NAs for graphing
graph_incidentmap$Total[graph_incidentmap$Total == 0] <- NA

## Make graph
ggplot() + 
  geom_map(data = graph_incidentmap, aes(map_id = FIPS, fill = Total), map = countymap) +
  scale_x_continuous( limits = c( -125 , -100 ) , expand = c( 0 , 0 ) ) +
  scale_y_continuous( limits = c( 30 , 50 ) , expand = c( 0 , 0 ) )  +
  labs(fill="Total Incidents") + coord_map() +
  scale_fill_gradient(low="gray80", high="gray28", na.value="white") +
  geom_path(aes(x=long, y=lat, group=group), data=cborders, color='gray14') + 
  geom_path(aes(x=long, y=lat, group=group), data=west_coast, color='black') + 
  theme_bw() + theme_void()


##################################################################
###                    Make map of sheriffs                    ###
##################################################################

## Use file "SheriffMap.dta"

## Merge in sheriff data
SheriffMap$FIPS <- as.numeric(as.character(SheriffMap$FIPS))
colnames(SheriffMap)[2] <- "sheriffdummy"
graph_sheriffmap <- merge(countymap,SheriffMap,by="FIPS", all = T)
graph_sheriffmap$sheriffdummy <- as.numeric(graph_sheriffmap$sheriffdummy)

## Make graph
ggplot() + 
  geom_map(data = graph_sheriffmap, aes(map_id = FIPS, fill = sheriffdummy), map = countymap) +
  scale_x_continuous( limits = c( -125 , -100 ) , expand = c( 0 , 0 ) ) +
  scale_y_continuous( limits = c( 30 , 50 ) , expand = c( 0 , 0 ) )  +
  labs(fill="") + coord_map() +
  scale_fill_gradient(low="gray80", high="gray28", na.value="white") +
  geom_path(aes(x=long, y=lat, group=group), data=cborders, color='gray14') + 
  geom_path(aes(x=long, y=lat, group=group), data=west_coast, color='black') + 
  theme_bw() + theme_void() +
  theme(legend.position="none")

#######################################################################
###            Make map of federal land ownership                   ###
#######################################################################

## Use the pct_federal.xls file for this
## Note for importing: the first row is data, not variable labels

pct_federal$pct_federal <- pct_federal$X__6
pct_federal$FIPS <- as.numeric(pct_federal$X__1)

## Merge in land ownership data
graph_land <- merge(countymap,pct_federal,by="FIPS", all = T)

## Need to have zeros become NAs for graphing
graph_land$pct_federal[graph_land$pct_federal == 0] <- NA

## Make new base map that is the whole country 
state_map_full <- map_data("state")
county_map_full <- map_data("county")

## Make graph
margin=margin(0, 7, 0, 7) 

ggplot() + 
  geom_map(data = graph_land, aes(map_id = FIPS, fill = pct_federal), map = countymap) +
  labs(fill="Percent Federal Land") + coord_map() +
  scale_fill_gradient(low="gray90", high="gray20", na.value="white") +
  geom_path(aes(x=long, y=lat, group=group), data=state_map_full, color='black') +
  theme_bw() + theme_void()

